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ABSTRACT 

We use three dimensional magnetohydrodynamic simulations, in a pseudo-Newtonian potential, to 
study geometrically thin accretion disc flows crossing rms, the marginally stable circular orbit around 
black holes. We concentrate on vertically unstratified and isothermal disk models, but also consider a 
model that includes stratification. In all cases, we find that the sonic point lies just inside rms, with 
a modest increase in the importance of magnetic field energy, relative to the thermal energy, observed 
inside the last stable orbit. The time-averaged gradient of the specific angular momentum of the flow, 
(dl/dr), is close to zero within r^g, despite the presence of large fluctuations and continuing magnetic 
stress in the plunging region. The result that the specific angular momentum is constant within is 
in general agreement with traditional disk models computed using a zero-torque boundary condition at 
the last stable orbit. 

Subject headings: accretion, accretion disks - black hole physics - MHD - stars: neutron - 
hydrodynamics - instabilities 



1. INTRODUCTION 

The existence of a marginally stable orbit is a distinc- 
tive feature of accretion flows onto black holes. Within 
rms, stable circular orbits do not exist, and gas inevitably 
plunges on a short timescale into the hole. The location of 
the marginally stable orbit, which lies at rms — 6GM/c^ 
for a non-rotating hole, varies strongly with the spin pa- 
rameter a for Kerr black holes. This variation, together 
with the assumption of a large change in the emission prop- 
erties of the gas interior and exterior to rms , is central to 
all attempts to measure a from observable quantities for 
both galactic (Zhang, Cui & Chen 1997) and supermassive 
black holes (Iwasawa et al. 1996; Dabrowski et al. 1997; 
Bromley, Chen & Miller 1997). Neutron stars, if they are 
sufficiently compact, could also possess a last stable orbit, 
though we will not consider this possibility further here. 

The presence of gas on plunging orbits within rms can 
have observable consequences. It potentially modifies 
the profile of the relativistic iron Ka line that has been 
observed in the X-ray spectra of some Seyfert galaxies 
(Tanaka et al. 1995; Reynolds & Begelman 1997; Young, 
Ross & Fabian 1998), and may create detectable absorp- 
tion signatures (Nandra et al. 1999). However, it has 
generally been believed that gas within rms does not have 
any dynamical mvpoTtance, since the infall rapidly becomes 
supersonic a short distance inside the marginally stable 
orbit. Whatever complexities may develop in the flow 
subsequently are then causally disconnected from the disk 
at larger radius. This line of reasoning suggests that for 
modelling the disk (by which we mean the region of the 
flow outside the last stable orbit), it suffices to impose a 



zero-torque boundary condition at rms and ignore the inte- 
rior region altogether. For the details of disk models con- 
structed in this manner, we refer the reader to Abramowicz 
& Kato (1989), and references therein. 

This simple hydrodynamic picture would be oversimpli- 
fied if magnetic fields inside rms were strong enough to 
maintain a connection between the plunging region and 
the disk. A recent analysis by Krolik (1999) showed that 
if magnetic fields, generated in the disk, remain frozen 
into the gas as it inspirals within rms, then the energy 
density in the fields can become comparable to the rest- 
mass energy density of the flow. Related ideas have also 
been proposed by Gammie (1999). In addition to altering 
the zero-torque boundary condition for the disk, the pres- 
ence of extremely strong fields in the plunging region could 
have several potentially observable consequences (Agol & 
Krolik 2000). For example, the radiative efficiency of the 
disk, e = L/Mc^, could be substantially increased. The 
analysis of Krohk (1999), however, depends upon a split 
of the accretion flow into two distinct regions, a disk at 
f > Tins in which magnetic dynamo processes (Balbus & 
Hawley 1991; Brandenburg et al. 1995; Stone et al. 1996; 
for a review see e.g. Hawley & Balbus 1999) maintain 
a turbulent state, and a plunging region at r < rms in 
which it is assumed that magnetic flux is frozen into the 
fluid. This split is clearly an approximation, because the 
same turbulent processes that act to generate and destroy 
magnetic flelds in the disk itself will continue to operate 
until some finite distance inside the last stable orbit. This 
may alter the conclusion that a rapid growth in the rela- 
tive importance of magnetic fields is inevitable (Paczynski 



^CITA, University of Toronto, McLennan Labs, 60 St George Street, Toronto, Ontario M5S 3H8, Canada 

^Present address: School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews, Fife, KY16 9SS, UK 
^Hubble Fellow 



2 



Accretion flows crossing the last stable orbit 



2000). The analytic analysis further assumes that the flow 
is axisymmetric and steady-state. This is also only ap- 
proximately the case in a turbulent accretion flow, and 
needs to be tested via numerical simulations. 

In this paper, we employ three dimensional magneto- 
hydrodynamic (MHD) simulations to study the proper- 
ties of the accretion flow as it crosses ry^s- A pseudo- 
Newtonian potential appropriate for a non-rotating black 
hole (Paczynski & Wiita 1980) is used to mimic the rela- 
tivistic effect of a last stable orbit within a non-relativistic 
MHD code. Within this approximation, simulating the 
evolution of magnetic flelds within the plunging region is, 
if anything, expected to be easier than following the devel- 
opment of turbulence in the disk at larger radii. Within 
the disk, orbits arc stable, and thus in the absence of tur- 
bulence a held loop in the plane of the disk will always 
be sheared out in azimuth until numerical reconnection 
sets in at the grid scale. Once into the plunging region 
there is only a finite (but large) amount of shear before 
the infalling gas reac;lies the black hole. 

The accretion flow within r^s has already been studied 
by Hawley (2000), who presented global MHD simulations 
of accretion tori with a geometry similar to that of pop- 
ular ADAF (Narayan & Yi 1994) and ADIOS (Blandford 
& Begelman 1999) models for radiatively inefficient fiows. 
This geometry is thought to be appropriate at relatively 
low accretion rates for both supermassive and stellar mass 
black hole systems (for observational support see, e.g., Gil- 
fanov, Churazov & Revnivtscv 1999). The specific model 
considered was geometrically fairly thick at large radius 
{r '> r^ns), although near the last stable orbit the im- 
portance of pressure forces, and the vertical scale height, 
was smaller. Indeed, at r = r^s, the relative scale height 
was {h/r) ~ 0.1, which is similar to the value used in our 
simulations. Significant magnetic stresses within rms were 
obtained in these simulations, and in subsequent higher 
resolution simulations of the same geometry (Hawley & 
Krolik 2000). 

Existing work suggests that the conditions at rms may 
differ for geometrically thin and thick disks (Popham & 
Gammie 1998). We therefore emphasize that in this paper 
we address only the case of geometrically thin disks. The 
existence of a thin disc implies that radiative cooling must 
be efficient, and this allows us to further simplify the prob- 
lem by adopting an isothermal equation of state in lieu of 
solving the energy equation. Thin disks are the expected 
mode of accretion at high M in both Active Galactic Nu- 
clei, and in stellar mass Galactic black hole sources. 

2. NUMERICAL METHODS 
2.1. The ZEUS hydrodynamics code 

We simulate the accretion flow using the ZEUS code 

(Stone & Norman 1992a,1992b; Clarke, Norman & Fielder 
1994; Norman 2000) to solve the equations of ideal mag- 
netohydrodynamics. ZEUS is an explicit Eulerian finite 
difference code, formally of second order accuracy, which 
uses an artificial viscosity to reproduce shocks. Advan- 
tages of the code for accretion disk applications include a 
fiexible choice of gridding and co-ordinate systems, and al- 
gorithms (Norman, Wilson & Barton 1980) that minimise 
spurious diffusion of angular momentum relative to mass. 
The analyis of the behavior of magnetic fields in the 



plunging region by Krolik (1999) is independent of the ex- 
istence of vertical stratification in the disk. For our initial 
simulations we therefore adopt the computationally easi- 
est option, and consider a vertically unstratified disk (i.e. 
one where the vertical component of gravity is artificially 
set to zero), and an isothermal fluid, where the presure P 
is given by, 

P = pcl (1) 

where Cg is the sound speed and p is the density. For our 
standard model we choose Cg such that the ratio of the 
sound speed to the circular orbital velocity at the last sta- 
ble orbit, Cs/v^ — 0.08. In an unstratified simulation, of 
course, Cg plays no role in setting the vertical scale height 
of the flow. However, the ratio Cg/v^, which in a real disk 
is related to the relative disk scale height via h/r ~ Cs/v^, 
still determines the importance of radial pressure forces for 
the disk structure. The simulated fiows are 'thin' in the 
sense that radial pressure forces, which scale as {h/r)^, are 
small compared to the gravitational force. A model with 
smaller sound speed was also investigated. For better com- 
parison with the simulations of Hawley (2000) , we also ran 
a variant including the effect of vertical stratification. 

Standard choices for the Courant number, Co = 0.5, 
and coefficient of artificial viscosity, C2 = 2.0, were used 
throughout. Second order advection was used for all quan- 
tities. The only code scaling which needs to be mentioned 
is that for time. We adopt units in which the orbital period 

^t Tins, ^ms — 7.T. 

2.2. Paczynski- Wiita potential 

ZEUS is a Newtonian fluid dynamics code which does 

not model any of the effects of special or general relativ- 
ity. Within this framework, we use a pseudo-Newtonian 
gravitational potential (Paczynski & Wiita 1980), 



where rg = 2GM/c^, to model what is expected to be the 
dominant relativistic effect around a non-rotating black 
hole - the existence of an innermost stable orbit at rms = 
QGM/c^. 

2.3. Initial and boundary conditions 

We simulate a wedge of the disk extending 30° in az- 
imuth in cylindrical polar geometry, (r, z, (j)). This limita- 
tion on the azimuthal extent of the domain has the unde- 
sirable effect of eliminating larger scale azimuthal modes of 
the magnetic field, which contribute significantly to the to- 
tal power in global disk simulations (e.g. Armitage 1998). 
The local properties of disk turbulence (for example the 
relative strengths of radial and azimuthal field compo- 
nents), however, which are probably of more importance 
for this problem, are less affected by the size of the az- 
imuthal domain. 

The boundary conditions are periodic in cj), and set to 
outflow at both rmin and rmax- Outflow boundary condi- 
tions arc implemented in ZEUS by setting zero gradients 
for all flow variables at the boundary. We note that this 
choice of boundary conditions allows a significant toroidal 
magnetic field to be advected across the inner boundary 
- we have not attempted to impose Newtonian analogs 
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of the general relativistic boundary conditions for mag- 
netic fields at the event horizon of the black hole. From a 
numerical standpoint, the use of outflow boundary condi- 
tions is desirable to ensure that the simulations do indeed 
model magnetic instabilities, rather than purely hydrody- 
namic instabilities to which relatively narrow annuli, with 
reflecting boundaries, are known to be susceptible (e.g. 
Narayan & Goodman 1989). We use a uniform grid in 
both (f) and z. In the radial direction, a uniform grid inte- 
rior to Tnis is matched smoothly onto a grid at larger radii 
for which Ar^+i = kAvi, with fc > 1 a constant. This con- 
centrates resolution in the region of greatest interest near, 
and within, the marginally stable orbit. 

The initial conditions for the calculation are intended 
as thin disk analogs of the toroidal configuration used by 
Hawley (2000), with the additional assumption of isother- 
mality. We take a gaussian surface density profile, cen- 
tered at r = 2ri„s) with an azimuthal velocity appropriate 
to Keplerian rotation. We relax this profile in a prelimi- 
nary non-magnetic one dimensional calculation, to ensure 
that it is an accurate numerical equilibrium state. The 
resulting profile of p and is then mapped into three 
dimensions, and a seed magnetic field added. To evalu- 
ate the sensitivity of the results to the boundary condi- 
tions, three unstratified simulations were run, with dif- 
ferent choices of seed field and associated boundary con- 
ditions in z. Two simulations were run with an initially 
vertical field imposed at all radii where the surface density 
exceeded a threshold value, taken to be approximately a 
quarter of the maximum value of the surface density. Pe- 
riodic boundary conditions for all variables were imposed 
in the z direction. For the first simulation (the 'standard' 
run) the radial and vertical boundaries were at rmin = 1, 
^■max = 5, and z = ±0.3, in code units where rms = 1-5- 
The initial field strength (in regions seeded with field) 
was such that the ratio of gas pressure to magnetic pres- 
sure, Pz = 5000. The numerical resolution was = 150, 
Hz = 32 and = 40 grid points, with 30 of the radial 
grid points interior to rms- The second run simulated a 
cooler disk (with a sound speed half the previous value), 
in a larger computational volume, bounded by 1 < r < 10 
and z — ±0.25. This run used Ur = 210, Uz = 32 and 
= 48 grid points, and an initial seed field /3z = 800. 
Finally, a simulation was run with parameters similar to 
the first, except starting with an initially azimuthal field, 
with (]^ — 100. For this simulation the vertical boundary 
conditions were chosen to be refiecting, with the vertical 
components of both the velocity and the magnetic field set 
to zero on the boundaries. The resolution for this run was 
120 X 32 X 40 grid points. 

Simulations that include vertical stratification are obvi- 
ously more realistic, but they are also more demanding of 
computational resources, because the development of mag- 
netically dominated low density regions at high \z\ places 
severe restrictions on the timestep. Numerical tricks can 
be used to mitigate this problem (Miller & Stone 2000), 
but for this paper we adopted the simpler approach of con- 
sidering a volume containing only a small number of disk 
scale heights. We therefore reran the standard model, in- 
cluding the vertical component of gravity, in a domain 
boimded hy z = ±0.5. This admits only a couple of scale 
heights at rms- As in the standard model, periodic bound- 
ary conditions were used for the magnetic field in the ver- 



tical direction. For this run we used n,. = 150, = 48 
and ricj, = 40 grid points. 

It is worth noting at the outset that these simulations 
are not, and are not intended to be, realizations of the 
same physical situation, and thus the results will differ 
between them. An initially vertical field gives the fastest 
possible growth rate of the Balbus-Hawley (1991) insta- 
bility. Additionally, a non-zero vertical flux boosts the 
strength of MHD turbulence and increases the effective 
Shakura-Sunyaev (1973) a parameter obtained (Hawley, 
Gammie & Balbus 1995). Substantially longer timescales 
are required to reach saturation with an initially azimuthal 
fleld. We have run both simulations that used the smaller 
computational domain until a significant fraction (around 
a half) of the initial mass had been accreted, and compare 
the results at this late stage when the fiow in the inner 
regions of the disk, near rms, has reached an approximate 
steady state. 

3. RESULTS 

Fig. 1 shows the growth in the magnetic energy of the ra- 
dial magnetic field component in the standard unstratified 
model with an initially vertical (z direction, f3z = 5000) 
magnetic field. Around 10 orbits of evolution are required 
to reach a nonlinear state in the inner regions of the flow, 
during which time the radial magnetic field energy grows 
exponentially. The growth rate is consistent with the ex- 
pected growth rate of the Balbus-Hawley instability in 
this field geometry (Balbus & Hawley 1991), evaluated at 
the smallest radius that has been seeded with magnetic 
field. Subsequent to the instabilities saturating, the field 
energy fluctuates but remains roughly constant, before fi- 
nally starting to decline as a significant fraction of the disk 
is accreted. At the point the run was stopped, a,t t = 200, 
65% of the initial mass had been accreted. The magnetic 
field was dominated by the toroidal component, with the 
ratio of magnetic field energy in the z, r and (j) fields being 
approximately 1:2:30. 

The simulation seeded with an azimuthal magnetic field 
shows similar behavior, but with much slower growth of 
the magnetic field energy, despite the initially higher ratio 
of the magnetic energy to the thermal energy. This config- 
uration was run up to t = 600, by which time just under 
40% of the initial disk mass had been accreted. 

Fig. 2 shows the evolution of the disk surface density 
with time, for the standard nm with an initial vertical 
field. The evolution differs somewhat from the standard 
diffusive evolution of a viscous annulus (e.g. Pringle 1981), 
due to both the varying amount of time required before 
the disk at different radii becomes fully turbulent, and 
due to the development of supersonic infall interior to 
Tms • However the general trend towards a broadening sur- 
face density profile, that remains relatively smooth, is re- 
covered. Spatial fiuctuations in E, shown in Fig. 3, are 
strongly sheared and thus predominantly azimuthal in ex- 
tent. When normalised to the mean surface density as a 
function of radius, as in the figure, there is no obvious 
visual indication of the location of the marginally stable 
orbit. 

The radial and azimuthal velocities as a function of ra- 
dius obtained in the standard run are shown in Fig. 4. 
A small radial velocity within the disk itself transitions 
smoothly into rapid infall within the marginally stable or- 
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bit. The sonic point hcs somewhat inside r^s (within a 
radial distance ^ h, the disk scale height that would cor- 
respond to the assumed sound speed). Note that the radial 
velocity slightly outside the marginally stable orbit is al- 
ready beginning to increase in magnitude, in advance of 
reaching the actual infall region. Beyond r « 2r^s, the ra- 
dial velocity is outward. This is a consequence of the initial 
surface density profile and the use of a free outer boundary 
condition. We have also experimented with runs in which 
mass was continually injected across the outer boundary, 
and the disk allowed to evolve until a quasi-steady state 
was achieved. These runs yielded similar results to those 
reported here - in particular there was no strong ampli- 
fication of magnetic field observed interior to rms- How- 
ever, with inflow boundary conditions at rout it is hard 
to be sure that purely hydrodynamic instabilities, which 
are known to occur within relatively narrow annuli when 
the boundary conditions are reflecting, do not contami- 
nate the results. We therefore concentrate attention on 
the current simulations which follow a ring of gas which is 
allowed both to accrete and to spread to larger radii. 

As expected for a thin disk, in which radial pressure 
forces are small compared to the gravitational force, the 
azimuthal velocity in the disk, v^, is very close to the Kc- 
plerian value for orbits in a Paczynski-Witta potential. 

The presence of magnetic fields in the flow potentially 
allows the fiow interior to the sonic point to communi- 
cate with the exterior disk. Gcnerically, except in unusual 
circumstances, MHD turbulence driven by the Balbus- 
Hawley instability saturates at a level where the mag- 
netic pressure in the disk is substantially less than the 
thermal pressure, typically by one or two orders of magni- 
tude (Hawley, Gammie & Balbus 1995; Brandenburg et al. 
1995; Stone et al. 1996). Stronger fields, relative to the 
gas pressure, are likely to occur in the disk corona (Miller 
& Stone 2000). As shown in Fig. 5 for our standard un- 
stratified simulation, the mean Alfven speed in the disk, 
VA = y/ {B'^/AiTp), remains less than the sound speed at 
all radii. This is even more true for the simulation with 
the initial (p field, not shown in the figure, since the sat- 
uration level of the magnetic fields in that run is signifi- 
cantly below the value obtained in the initial vertical field 
runs. Similarly, for the vertically stratified version of the 
standard run, the mean Alfven speed is smaller than the 
sound speed at all radii, and is of similar magnitude to 
the unstratified case. For this run, the magnetic fields and 
the Alfven speed near the disk midplane, also plotted in 
Fig. 5, are significantly greater than the mean value in the 
plunging region, though not by a large factor. 

These results do not, however, exclude the possibility 
that the plunging region can exert a torque on the disk 
at the last stable orbit. In the presence of turbulence, 
some regions of the fiow have relatively stronger mag- 
netic fields (or lower density), with correspondingly higher 
Alfven speed. As a result, the critical point can penetrate 
further into the plunging region than one would estimate 
based solely on the mean flow properties. As shown in 
Fig. 5, for both the stratified and the unstratifed runs, we 
find that the peak Alfven speed at any radius is substan- 
tially greater than the mean value. Although it does not 
rise as rapidly as the absolute value of the radial veloc- 
ity, it is therefore possible, at least in principle, for the 



plunging region some distance inside the last stable orbit 
to be in communication with the disk at larger radii. In 
our simulations this distance is still only ^ h, but whether 
this result also holds for cooler flows is currently unknown. 

Fig. 6 shows the ratio of the magnetic energy to the 
thermal energy for the standard run, and for the run with 
an initially azimuthal seed magnetic field, as a function of 
radius. A single timeslice from the simulation produces a 
rather noisy estimate of this quantity, so we plot the av- 
erage from 5 independent epochs close to the end of the 
run {t = 150 to t = 200 for the z field run, and t = 440 
to t = 600 for the field case). The ratio (B'^/Snpcl) is 
around 0.05 - 0.1 for the initial z field nm with (3^ — 5000, 
and is modestly greater for the run using a larger compu- 
tational domain and an initial Pz = 800. These vahies are 
consistent with the values (0.1 - 0.2) obtained by Hawley, 
Gammie & Balbus (1995) from unstratified local simula- 
tions in the shearing-box geometry, at comparable initial 
Pz = 3200. As expected, the fields are significantly weaker 
in the simulation with an initial 4> field and non-periodic 
boundary conditions in z. In this case we obtain a ratio 
of magnetic to thermal energy in the 1 — 2 x 10~^ range. 
In all the runs, modest field amplification interior to the 
marginally stable orbit is observed. The fields remain well 
below equipartition with the thermal energy, while the ki- 
netic energy is of course orders of magnitude greater still. 

Angular momentum transport in the disk is dominated 
by magnetic, rather than fluid, stresses. We also plot in 
Fig. 6 a measure of the magnetic stress, normalised to the 
gas pressure, 



which in the disk is just the magnetic contribution to the 
Shakura-Sunyaev a parameter. Typical values obtained 
are a few x 10~^ for the standard simulation, and some- 
what less than 10~^ for the simulation with an initial mag- 
netic field in the (p direction. From the figure, it can also be 
seen that there is continuing magnetic stress, and a rise in 
Q^mag) within the plunging region, where the density drops 
rapidly. In these plots there is evidently nothing special 
about the location r = rms, and the magnetic torque does 
not vanish there. The existence of this stress inside rms, 
which is one of the predictions of Krolik (1999), has also 
been seen in previous numerical simulations (Hawley 2000; 
Hawley & Krolik 2000). 

For the vertically stratified run, we have also checked for 
the presence of significant variations of Q!mag, and other 
quantities, with z. Interior to rms, the strength of the 
magnetic stress, normalised to the gas pressure, is found 
to be greater near the midplane of the disk than at larger 
z, but only modestly so. 

To test the eff'ect of this ongoing stress on the dynamics 
of the flow inside the last stable orbit, we plot in Fig. 7 
the specific angular momentum I of the flow as a function 
of radius. To reduce fluctuations in we average the spe- 
cific angular momentum over several timeslices taken from 
near the end of the simulations. For all the runs, /(r) in 
the disk (r > rms) is close to the Keplerian value for circu- 
lar orbits in the pseudo-Newtonian potential, apart from 
a deviation close to rms where the radial density gradient 
is becoming large. At a zero-torque boundary, we expect 
that dl/dr vanishes. There is some variation in I within 
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the marginally stable orbit, even in the averaged profiles, 
and weak indications of a small systematic decline towards 
smaller radii. This is consistent with the aforementioned 
presence of magnetic stress inside rms- However, it is clear 
from Fig. 7 that for both the unstratified and stratified 
simulations the specific angular momentum is close to flat 
within Tms- The magnetic stress that exists within ry^s 
does not appear to be sufficient to change I significantly 
within the last stable orbit. In this (restricted) sense, the 
numerical results are broadly consistent with the standard, 
purely hydrodynamic, picture of thin disk accretion accre- 
tion onto black holes (Abramowicz & Kato 1989), which 
assumes a zero-torque boundary condition at the last sta- 
ble orbit. 

The highest resolution of the runs discussed here is ob- 
tained in the vertically stratified simulation. For this run, 
we plot in Fig. 8 the specific angular momentum from 5 
timeslices, evenly spaced between t = 150 and t = 200. 
Significant fluctuations in I at any given radius are present. 
The influence of such fluctuations on the average structure 
of the disk outside Tms is unclear, though we speculate that 
they could lead to significant changes in, for example, the 
radiative efficiency of the flow. However, in none of the 
slices do we see evidence for the clear and continuing de- 
cline in I, interior to the marginally stable orbit, that was 
obtained in the simulations of Hawley (2000). 

4. DISCUSSION 

In this paper, we have used MHD simulations to study 
the transition between a geometrically thin accretion disk, 
in which inflow is driven by the rate at which turbulence 
can transport angular momentum outwards, and the un- 
stable plunging region interior to the marginally stable 
orbit. We find that in many respects the transition re- 
sembles that expected on the basis of previous analytic 
and one-dimensional numerical calculations (Paczynski & 
Bisnovatyi-Kogan 1981; Muchotrzeb & Paczynski 1982; 
Muchotrzeb 1983; Matsumoto et. al. 1984; Abramow- 
icz & Kato 1989). In particular, we find no evidence, at 
least in pseudo-Newtonian simulations at the current reso- 
lution, for the extremely strong growth in the importance 
of magnetic fields (relative to the thermal or rest-mass 
energy of the flow) and associated dynamical effects dis- 
cussed by Krolik (1999). Indeed, in these thin disk sim- 
ulations, the average specific angular momentum is close 
to flat at and within the last stable orbit. In this respect, 
the numerical results, at least for the time averaged state, 
appear to be consistent with disk models computed using 
a zero-torque boundary condition at the last stable orbit. 
This differs from the results of some previous simulations 
(Hawley 2000; Hawley & Krolik 2000), in which an unmis- 
takeable decline in the specific angular momentum inside 
Tins was obtained. Although the conditions at rms are ex- 
pected to vary between thin and thick disks, the model 
computed by Hawley (2000) and Hawley & Krolik (2000) 



is thin enough at r^s that l{r) is close to Keplcrian. We 
speculate that differences in the equation of state, spatial 
domain, or numerical resolution could be responsible for 
the discrepancy. Further simulations investigating these 
possibilities are in progress. 

Although we have not found any strong dynamical ef- 
fects associated with magnetic fields interior to the last 
stable orbit, we do see evidence for the essential ingre- 
dient - ongoing magnetic stresses in the plunging region 
- of recent models that have questioned the validity of a 
zero-torque boundary condition for black hole accretion 
(Gammie 1999; Krolik 1999; Agol & Krolik 2000). The 
differences between these simulations, and those presented 
by Hawley (2000) and Hawley & Krolik (2000), are thus 
quantitative rather than qualititive in nature. We woiild 
therefore emphasize that further improvements in the sim- 
ulations are required to investigate the behavior of the ver- 
tically stratified flows more robustly, to extend the simu- 
lations to test the evolution of cooler flows crossing the 
last stable orbit, and to explore the continuum between 
geometrically thin and thick disks. 

The simulations reported here are non-relativistic, and 
cannot - even crudely - be extended to model any of the 
additional complexity that arises if the accretion flow is 
onto a Kerr blac;k hole. Considerable progress has been 
made in the development of numerical methods for gen- 
eral relativistic hydrodynamics and magnetohydrodynam- 
ies (e.g. Font 2000), but these methods are probably not 
yet able to follow a turbulent disk flow transiting the last 
stable orbit. We note, however, that the general relativis- 
tic MHD calculations that have already been reported do 
show important differences with non-relativistic analogs 
(Koide, Shibata & Kudoh 1999; Meier 1999; Koide et al. 
2000), and this should be borne in mind as an impor- 
tant caveat to the results presented here. The dynam- 
ics of the inner disk could also be modified by magnetic 
fields exterior to the disk, either in a disk corona (Miller 
& Stone 2000; Hawley & Krolik 2000), or by fields finking 
the disk to a spinning black hole (Livio 1999; Blandford 
1999). Given the diverse range of variability observed in 
accreting black hole systems, further exploration of such 
ideas is clearly warranted. 
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Fig. 1. — Energy in the radial component of the magnetic field, integrated over the simulation volume, as a function of time. The units on 
the vertical axis are arbitrary. The initial magnetic field for this simulation was in the z direction, and the vertical boundary conditions were 
periodic. The dashed line shows the expected growth rate of the radial magnetic field energy density, for the most unstable (largest f2) mode. 
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Fig. 2. — Evolution of the disk surface density profile in the simulation with an initially vertical magnetic field geometry. Prom top down, 
the slices are plotted at t = 0, t = 50, t = 100, t = 150, and t = 200. More than half of the mass has been accreted by the end of the 
simulation. 
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Fig. 4. — Radial and azimuthal velocity at f = 200 from the standard simulation. The short dashed curve in the lower panel shows the 
Keplerian velocity of circular orbits in the pseudo-Newtonian potential, in units where c = 1. 
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Fig. 5. — Wave propagation speed as a function of radius, evaluated at t = 200 from the standard simulation with and without vertical 
stratification. In the upper panel, for the unstratified run, the horizontal dotted line shows the sound speed, the short dashed line the mean 
Alfven speed as a function of radius. The negative of the radial velocity is plotted as the solid curve. The lower panel shows the same 
quantities plotted for the stratified simulation, with the long dashed line showing the mean Alfven speed in the flow near the disk midplane 
(|z| < 0.1). We also plot the peak Alfven speed at each radius (dot-dashed line), which exceeds the mean value by a substantial factor. 
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Fig. 6. — The upper panel shows the ratio of the magnetic energy to the thermal energy, as a function of radius, for the standard simulation 
with an initial z field (dashed line), and for the run with an initial (p field (solid line). In both cases, several timeslices from near the end of 
the runs have been averaged together to reduce the magnitude of the fluctuations. The lower panel shows the effective Shakura-Sunyaev a 
parameter derived from the magnetic torques. 
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Fig. 7. — The specific angular momentum I as a function of radius. Tfie solid line shows the result for the unstratified simulation with 
an initial <f) magnetic field, the short dashed and long dashed curves for unstratified simulations with an initial z field but different sound 
speeds and seed field strength. The dot-dashed line shows the result for the vertically stratified simulation, evaluated in the disk midplane. 
Note that the simulation with reduced sound speed (long dashed line) evolves more slowly, so that the averaging of several timeslices has not 
eliminated substantial fluctuations present at larger radius outside f'ms. For all these models, dl/dr is close to zero at rms- The dotted curve 
shows the specific angular momentum of circular orbits in the pseudo-Newtonian potential. 
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FlH. 8. The specific angular momentum / as a function of radius, plotted from five timeslices of the stratified run. As in Fig. 7, / is 
evaluated in the disk midplane. The dotted curve shows the specific angular momentum of circular orbits in the pseudo-Newtonian potential. 



